Black-box superconducting circuit quantization 
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We present a semi-classical method for determining the effective low-energy quantum Hamiltonian 
of weakly anharmonic superconducting circuits containing mesoscopic Josephson junctions coupled 
to electromagnetic environments made of an arbitrary combination of distributed and lumped ele- 
ments. A convenient basis, capturing the multi-mode physics, is given by the quantized eigenmodes 
of the linearized circuit and is fully determined by a classical linear response function. The method 
is used to calculate numerically the low-energy spectrum of a 3D-transmon system, and quantitative 
agreement with measurements is found. 
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Superconducting electronic circuits containing nonlin- 
ear elements such as Josephson junctions (JJs) are of 
interest for quantum information processing [1, 2], due 
to their nonlinearity and weak intrinsic dissipation. The 
discrete low-energy spectrum of such circuits can now 
be measured to a precision of better than one part per 
million [.'-!]. The question thus naturally arises of how 
well one can theoretically model such man-made artifi- 
cial atoms. Indeed, increasing evidence indicates that 
due to increased coupling strengths [4], current models 
are reaching their limits [5-9] and in order to further our 
ability to design, optimize and manipulate these systems, 
developing models beyond these limits becomes neces- 
sary. This is the goal of the present work. 

An isolated ideal JJ has only one collective degree of 
freedom: the order parameter phase difference Lp across 
the junction. The zero-temperature, sub-gap physics of 
this system, with Josephson energy Ej and charging en- 
ergy E'c is described by the Cooper-pair box Hamilto- 
nian 
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where N is the Cooper-pair number operator conjugate 
to if> and Ng an offset charge. This model is exactly solv- 
able in terms of Mathieu functions [10, 11]. The crucial 
feature that emerges from this solution is that the charge 
dispersion, i.e. the maximal variation of the eigenenergies 
with Ng, is exponentially suppressed with Ej/Ec while 
the relative anharmonicity decreases only algebraically 
with a slow power-law in Ej/Ec- As a consequence, 
there exists a regime with Ej 3> Ec - the transmon 
regime - where the anharmonicity is much larger than the 
linewidth (e.g. due to fluctuation of the offset charge Ng), 
thus satisfying the operability condition of a qubit [12]. 
This is the regime of interest here. 

In order to be useful for quantum information pro- 
cessing tasks, several Josephson qubits must be made 
to controUably interact with each other and spurious in- 
teractions with uncontrolled (environmental) degrees of 
freedom must be minimized. In circuit quantum elec- 



~Eare cavi;t:y \ 
^ Dressed cavil;y 




^Dipole ^ntenna ' uyjo^epl^son /unction 




\ ' ^ 



FIG. 1. (Color online) Cartoon of a JJ at the center of a 
broadband dipole antenna inside a 3D microwave cavity. The 
presence of the antenna alters the geometry of the cavity- 
mode (full (red) curve) and a precise description requires the 
inclusion of many bare modes (dashed curves). 



trodynamics [2, 11, 13] (cQED), this is achieved by cou- 
pling the JJs to a common microwave environment with 
a desired discrete mode structure. So far such systems 
have mostly been described theoretically by models well 
known from quantum optics such as the single-mode 
Jaynes-Cummings model and extensions thereof [14]. 

When applied to superconducting circuits with multi- 
level artificial atoms, multi-mode cavities and increased 
coupling strengths [4, (>, 7] however, several technical and 
practical difficulties with these approaches arise. For ex- 
ample, capturing important effects of non-computational 
qubit states requires going to high orders in perturbation 
theory [ I 'l]. Also, determining the bare Hamiltonian pa- 
rameters, in terms of which these models are defined, is 
cumbersome and requires iterating between experiment 
and theory. Perhaps even more important arc the short- 
comings of the traditional approaches in dealing with 
the multiple modes of the cavity. Indeed high-energy, 
off-resonant cavity modes have already been measured 
to contribute substantially to the inter-qubit interaction 
strength [8, 15] and, via the multi-mode Purcell effect. 
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also to affect the coherence properties (relaxation and 
dcphasing) of the qubits [5]. Attempts at including this 
multi-mode physics in the standard models however, lead 
to difficulties with diverging series and QED renormal- 
ization issues [8, 10], which to the best of our knowledge 
remain unresolved. Fig. 1 illustrates the origin of the 
problem with the example of a JJ inside a 3D cavity (3D- 
transmon) [.3] . The presence of a relatively large metallic 
dipole antenna [17] can strongly alter the geometry of 
the cavity modes. This essentially classical effect, can be 
accounted for precisely only by including a sufficiently 
large number of bare modes. 

In contrast, we propose to start by considering the cou- 
pled but linearized problem in order to find a basis that 
incorporates the main effects of the coupling between 
multi-level qubits and a multi-mode cavity and then ac- 
count for the weak anharmonicity of the Joscphson po- 
tential perturbatively. The crucial assumption made here 
is that charge dispersion effects can be safely neglected. 
This is reasonable given that in state-of-the art imple- 
mentations of transmon qubits [3, 18], charge dispersion 
only contributes a negligible amount to the measured 
linewidths. Previous work discussed the nonlinear dy- 
namics of a JJ embedded in an external circuit classi- 
cally [19]. Here we go one step further and show how the 
knowledge of a classical, in principle measurable, linear 
response function lets us quantize the circuit, treating 
qubits and cavity on equal footing. 

Single junction case. We consider a system with a JJ 
with bare Joscphson energy Ej and charging energy Ec, 
in parallel with a linear but otherwise arbitrary electro- 
magnetic environment as depicted in Fig. 2 (a). Neglect- 
ing dissipation, the unbiased junction alone is described 
by the Hamiltonian (1). At low energies, when Ej ^ Ec, 
quantum fluctuations of the phase ip across the junction 
are small compared with tt and, as emphasized in the in- 
troduction, the probability of quantum tunneling of the 
phase between minima of the cosine potential is negligi- 
bly small. It is then reasonable to expand the latter in 
powers of ip, thus obtaining the approximate circuit rep- 
resentation of Fig. 2 (b), in which the spider symbol [19] 



represents the purely nonlinear part and Lj 



'/Ej 



and Cj = e^/{2Ec) the linear parts of the Joscphson 
element. Here 4>o = H/{2e) is the reduced flux quantum. 
To leading order, the energy of the spider element is given 
by = -0oVV(24ij). 

A quantity of central importance in the following is 
the impedance Z{uj) of the linear part of the circuit de- 
picted in Fig. 2 (c). The latter is a complex meromorphic 
function and by virtue of Foster's theorem [20, 21] can 
be synthesized by the equivalent circuit of parallel LCR 
oscillators in series shown in Fig. 2 (d). Explicitly 



M 



p=l 



1 



1 



(2) 




(d) 



Ri 



L2 
R2 



Rm 



Ci C2 Cm 

FIG. 2. (Color online) (a) Schematics of a JJ ((red) boxed 
cross) coupled to an arbitrary linear circuit (striped disk), 
(b) The Josephson element is replaced by a parallel combi- 
nation of: a linear inductance Lj, a linear capacitance Cj 
and a purely nonlinear element with energy Ej{l — cos(</5)) — 
{Ej /2)if'^ , represented by the spider symbol, (c) The linear 
part of the circuit shown in (b) is lumped into an impedance 
Z{u}) seen by the nonlinear element, (d) Foster-equivalent 
circuit (pole-decomposition) of the impedance Z{oj). 



where M is the number of modes [22] and we have 
adopted the electrical engineering convention of writing 
the imaginary unit as j = —i. This equivalent circuit 
mapping corresponds, in electrical engineering language, 
to diagonalizing the linearized system of coupled har- 
monic oscillators. The resonance frequencies of the lin- 
ear circuit are determined by the real parts of the poles 
of Z or more conveniently by the real parts of the ze- 
ros of the admittance defined as Y{uj) = Z{uj)~^, and 
for weak dissipation, i.e. Rp ^ Lp/Cp, are given 
by ujp = (LpCp)^^. The imaginary parts of the roots 
{2RpCp)~'^ , give the resonances a finite width. The ef- 
fective resistances are given by Rp = l/ReY{ujp) and 
the effective capacitances are determined by the fre- 
quency derivative on resonance of the admittance as 
Cp = (l/2)Imy (wp). Here and in the following the prime 
stands for the derivative with respect to frequency. Note 
that ImY''(a;) > [20]. Together this yields a compact 
expression for the quality factor of mode p: 



ujp lmY'{ujp) 
~2 Rer(wp) ■ 



(3) 



When applied to the mode representing the qubit, Eq. (3) 
gives an estimate for the Purcell limit on the qubit life- 
time T\ — Qqb/i^qb due to photons leaking out of the 
cavity. 

In order to derive the effective low-energy quantum 
Hamiltonian of the circuit, we next neglect dissipation 
{Rp — )> 00) and introduce the normal (flux) coordinates 
(\,p{t) = /pe^"f* + (/p)*e"-'''^''* associated with each LC 
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oscillator in the equivalent circuit. We can then imme- 
diately write the classical Hamiltonian function of the 
equivalent circuit as Hq = '^'^p'LiifpYi-^p)^^ fp^ where 
the subscript indicates that we consider the linear 
part of the circuit. Kirchhoff's voltage law implies that 
up to an arbitrary constant, = X^pli '?^p('')i where 
(j){t) — V{T)dT is the flux coordinate of the junction 
with voltage V{t). Note that by the second Joscphson 
relation, the order parameter phase difference is related 
to the latter via (p{t) = (j){t)/(po (modulo 2tt). 

Quantization is achieved in the canonical way [23, 24] 
by replacing the flux amplitudes of the equivalent oscil- 
lators by operators as 

with the dimensionless bosonic annihilation (creation) 
operators Op (aj,). Direct substitution yields the Hamil- 
tonian Hq = h^jpajjUp of M uncoupled harmonic 
oscillators (omitting the zero-point energies) and the 
Schrodinger operator of flux across the junction is 

M rr 

p=i 

We emphasize that the harmonic modes Op represent col- 
lective excitations of the linear circuit and their frequen- 
cies Wp are the equivalent of dressed oscillator frequencies. 
The coupling in the linear circuit is treated exactly and 
in particular no rotating wave approximation is used. 

The Hamiltonian of the circuit including the JJ is 
then H = Ho + ffni, where H^i = -(0)V(240o^ij) + 
O{{(l)/4>o)^). Physical insight may be gained by treating 
the nonlinear terms as a perturbation on top of Hq as- 
suming the eigenstates \ni, 712, . . . , of the latter with 

energies -£^ni^,n2, nj^j 

= J2i ^ihwi, to be non-degenerate. 
Considering only the leading order cj)'^ nonlinearity, one 
then obtains the reduced Hamiltonian 

H4 = H(^ + ^Y1 ^PP' ^P^P' ■ 
pp' 

Here fip = ajop and Hq = Hq + J^p ^p^p in- 
cludes a correction to the Lamb-shift given by Ap = 

- £■ - m^fl"^) ■ We have further in- 

troduced the generalized ^-shift Xpp' between modes p 
and p' . Clearly, ap = Xpp is the anharmonicity of the 
first excited state (self-Kerr) of mode p while Xpp' = Xp'p 
with p ^ p' IS the state-dependent frequency shift per 
excitation (cross-Kcrr) of mode p due to the presence of 
a single excitation in mode p' . Explicitly we find 

L C J 

Xpp ^ J y^-Ec, Xpp' ^ ^'^\/XppXp'p' ■ C^) 



Note that all modes acquire some anharmonicity due 
to the presence of the nonlinear J J. There is thus no 
strict separation of qubit and cavity anymore. Collo- 
quially, a mode with strong (weak) anharmonicity will 
be called qubit-like (cavity- like). Interestingly, in this 
lowest order approximation, the anharmonicity of mode 
p is seen to be proportional to the inductive partici- 
pation ratios [19] ip = Lp/Lj and inversely propor- 
tional to the capacitive participation ratio Cp = CpjCj. 
In the absence of a galvanic short of the junction in 
the resonator circuit, as is the case e.g. for a trans- 
mon qubit capacitively coupled to a cavity, it follows 
from the sum rule lim^j^o [^(^)/(i^)] = J2p^p = 
that «p < 1. Similarly, in the absence of any capaci- 
tance in series with Cj, it follows that Cp > 1, because 
lim^^o [jujZiLj)] = Cp-i = where Ce = Cj+Cy 
and C|| is the total capacitance in parallel with Cj. Hence 
we see that in this experimentally relevant case, the ef- 
fective anharmonicity of the qubit-like mode is always 
reduced as compared with the anharmonicity of the bare 
qubit given by — i?c [II]- Remarkably, in this approx- 
imation we find (see Eq. (7)) that the cross-Kerr shift 
between two modes is twice the geometric mean of the 
anharmonicities of the two modes. 

We emphasize that the above expressions do not how- 
ever account for higher order effects in anharmonicity 
such as the change of sign of the cross-Kerr shift ob- 
served in the straddling regime [If, 25]. Such effects are 
however fully captured by the full model H = Hq + H^i, 
which can be solved numerically. Remarkably, because 
the dressed modes already resum all the bare harmonic 
modes, typically only a few dressed modes M* ^ M 
need to be included for good convergence, thus consid- 
erably reducing the size of the effective Hilbert space, 
which scales as npii(-^p + 1) where Np is the maximal 
allowed number of excitations in mode p (e.g. Np = 1 in 
a two- level approximation). 

Charge dispersion. By assumption charge dispersion 
effects are neglected in the above approach. One may 
however ask how the charge dispersion of an isolated JJ is 
affected when the latter is coupled to a cavity. As in the 
Caldeira-Leggett model [2G], the coupling between the 
JJ and Harmonic oscillators suppresses the probability 
of flux tunneling and hence reduces charge dispersion of 
the qubit further. A simple estimate of the suppression 
factor is provided by the probability Pq of leaving the 
circuit in the ground state after a flux tunneling event 
and is found to be given by the "Lamb-Mossbauer" factor 

Pq w e~2 5:p^qb(2C^)/(''"p)^ where the sum excludes the 
qubit mode and Sq = Cj(})q/t is the charge (momentum) 
kick generated by a (j)Q flux slip through the J J of duration 
r and Cp = (l/2)Iniy (wp). Thus our assumption of 
neglecting charge dispersion of the qubit is well justified. 

Interestingly though, each eigcnmode of the system in- 
herits some charge dispersion. This effect, essentially a 
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i^c (GH 




1^02 


(GHz) 


aqb (MHz) 


X (MHz) 


Lj (nH) 


Cj (ff) 


7.77 


(7.763) 


8.102 


(8.105) 


15.33 


(15.333) 


-210 (-193) 


-90 


(-80.6) 


5.83 


7.6 


7.544 


(7 54') 


8.126 


(8.05) 


14.808 


(^A 830^ 


-280 f-24Q^ 


-30 


(-33.0) 


6.12 


9.2 


7.376 


(7.376) 


7.858 


(7.864) 


14.489 


(14.495) 


-264 (-257) 


-37.5 


(-38.7) 


6.67 


4.0 


7.058 


(7.045) 


8.005 


(8.023) 


13.788 


(13.794) 


-328 (-295) 


-13.2 


(-13.3) 


7.45 


5.2 


6.808 


(6.793) 


8.019 


(8.017) 


13.286 


(13.294) 


-330 (-293) 


-8 


(-8.4) 


7.71 


7.8 


6.384 


(6.386) 


7.832 


(7.823) 


12.45 


(12.449) 


-318 (-324) 


-5.4 


(-7.6) 


9.40 


0.34 



TABLE I. Low-energy spectrum (i^oi, I'c, V02), qubit anharmonicity (dqb) and state-dependent cavity shift (x) of six 3D- 
transmons. Results are shown in the format: experiment (theory). The theory values are obtained from a least square fit 
in Cj of the numerically computed lowest three energy levels of the 4>^ model. The fitted values of Cj are given in tlie last 
column. Their order of magnitude (a few femto-farads) agrees with estimates based on the sizes of the junctions. The Josephson 
inductances Lj are obtained from room-temperature resistance measurements of the junctions. 



consequence of hybridization, is of particular importance 
for applications such as quantum information storage in 
high-Q cavities coupled to JJs and is the subject of work 
in progress. 

Generalization to N junctions. The approach can be 
extended to circuits with multiple JJs connected in paral- 
lel to a common linear circuit. Details about the deriva- 
tion are given in the supplementary material [27] and we 
here only state the results. For N qubits, the resonance 
frequencies of the linear part of the circuit are deter- 
mined by the zeros of the admittance Yk{uj) = Zkk{'-^)^^ 
for any choice of reference port fc = 1, . . . , iV, where Z is 
the N X N impedance matrix of the linear part of the cir- 
cuit with a port being associated with each junction. The 
flux operators of the N junctions, with reference port fc, 
are given by = 1 , . . . , ) 




where Z^^ = 2/[ujplmYl{ojp)]. Note that the resonance 
frequencies are independent of the choice of reference 
port, while the eigenmodes do depend on it. In lowest 
order of PT and in the 0'* approximation, we find 

as well as the correction to the Lamb-shift Ap = G^pppp — 

^'^'^q l^qqPP- Here jiqq'ppi = X]s=l ^^^^T^Csg^sq'^sp'^sp' , 

and choosing the first port as the reference port (fc = 1), 
^"P ~ gill"'') \J '^ip • l^otice that the Cauchy-Schwarz 
inequality implies that |xgp| < 2y/a^7I^. Also, if q and 
q' refer to two different qubit-like modes, then Xqq' is 
a measure for the total interaction strength (cavity me- 
diated and direct dipole-dipolc coupling) between these 
two qubits. 

Comparison with experiment. As a demonstration of 
this method, we apply it to the case illustrated in Fig. 1 
of a single J J coupled to a 3D cavity [3]. The ad- 
mittance at the junction port y is a parallel combi- 
nation of the linearized qubit admittance and the ad- 
mittance Yc of the cavity-antenna system, i.e. ^(a;) = 



jujCj — i I [ujLj) + Yc(w). The junction is assumed to be 
dissipationless corresponding to a Purcell-limited qubit 
and ohmic losses of the cavity are included in Yc , which is 
complex. The Josephson inductance Lj is deduced from 
the measured junction resistance at room-tcmperature 
Rt, extrapolating it down to the operating tempera- 
ture [28] of 15 mk and using the Ambegaokar-Baratoff 
relation, Ej = h/S. / {^e^ Rt) . Cj - the only free parame- 
ter - is obtained by fitting the lowest three energy levels 
of the numerical solution of the 0^ model to the measured 
spectrum [ !]. Although Yc may in principle be obtained 
from current-voltage measurements, this is not practical 
in this system, where the antenna is hard to access non- 
invasively, being inside a closed high-Q cavity. Instead 
we use a finite element High Frequency Simulation Soft- 
ware (HFSS) and obtain Yc{l^) by solving the Maxwell 
equations numerically. Details on this simulation step 
are provided in the supplementary material [27]. 

From the zeros of the imaginary part of the admittance 
and their slopes we build and diagonalizc the (j)^ Hamilto- 
nian in a truncated Hilbcrt space, keeping in total three 
dressed modes (one qubit and two cavity modes) and 
allowing for maximally ten excitations per mode. The 
results of fitting the low-energy spectrum of six different 
samples are presented in Table I, where we also compare 
the predicted and measured qubit anharmonicities and 
shifts. We find agreement with the measured spectrum 
at the sub-per cent level and to within ten per cent with 
the measured anharmonicities and x-shifts. 

Conclusion and outlook. We have presented a sim- 
ple method to determine the effective low-energy Hamil- 
tonian of a wide class of superconducting circuits con- 
taining lumped or distributed elements. This method 
is suitable for weakly nonlinear circuits, for which the 
normal modes of the linearized classical circuit provide 
a good basis in the quantum case. For an N qubit 
system it requires only the knowledge of an A^ x 
(classical) impedance matrix. By working in a basis of 
dressed states, the parameters that appear in the Hamil- 
tonian incorporate much of the renormalization induced 
by the coupling between a multi-level artificial atom 
and a multi-mode resonator. Consequently, the number 
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of free parameters is considerably reduced as compared 
with standard models based on the Jaynes-Cummings 
paradigm expressed in terms of the experimentally in- 
accessible bare parameters. We have demonstrated the 
usefulness of this method in designing superconducting 
quantum information processing units by computing the 
low-energy spectrum of a 3D-transmon. Finally, this 
model may represent a suitable starting point for fu- 
ture investigations of the emerging ultra-strong coupling 
regime of cQED. 
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These notes provide further details on the HFSS simulation of the cavity admit- 
tance used to build the effective low-energy Hamiltonian in the black-box quantization 
approach to compare with the single junction experiment and on the black-box quan- 
tization method for the multi-qubit case. 



HFSS MODELING OF A 3D-TRANSMON 



As diseussed in the main text, the information about the spectrum of the quantum circuit, is 
encoded in the admittance at the port of the Josephson junction Y{u!) = Z(a;)~^. More precisely, 
it is sufficient to know the real roots and the derivative of Y at these points. 

Assuming that the size of the junction is negligibly small compared with the wavelength of 
the lower modes of the electromagnetic field in the cavity, it is appropriate to approximate the 
admittance of the linear part of the junction by a simple lumped element parallel LC oscillator 
with inductance Lj and capacitance Cj in parallel with the rest of the linear resonator. Hence the 
admittance can be decomposed as 

yH=j^Cj-^ + y,H, (lo) 

where Yc{uj) is the admittance of the system without the junction. The latter quantity can in 
principle be directly measured but in this particular design a measurement is not practical. Instead 
we simulate the classical system without the junction by solving Maxwell's equations numerically 
using HFSS. Fig. 3 shows a graphical representation of the different meshes used to represent 
the different elements of the cavity and antenna system. The smaller the element, the finer the 
mesh needs to be for accuracy and convergence. In this finite element simulation all metallic parts 
(Antenna and cavity boundaries made of pure aluminum), are treated as perfect conductors with 
zero resistance. In doing so, we neglect the kinetic inductance of the antenna and cavity. The 
finite London penetration depth of roughly A « 15 nm would lead to an effective increase of the 
cavity size and hence a decrease of the cavity frequency of about 10 kHz. Furthermore, the kinetic 
inductance of the antenna and wire connecting the two antenna pads to the Josephson junction 
can be estimated as 



2tanh(^) 



L I 



(11) 



where d w 100 nm is the thickness of the aluminum layer, L sa 1 mm is the total length and 
W sa 250 fim the width of the antenna and I ~ 34 /xm is the length and w ^ 1 /xm the width 
of the wire. With these numbers we obtain Lk ~ 1.6 • 10~^ nH, which is about three orders of 
magnitude smaller than (the linear part of) the Josephson inductance. A simple estimate shows 
that this would lead to a negative shift of the qubit resonance of only a few hundreds of kHz. 
These corrections are negligible at the current level of accuracy but can be easily included in the 
numerical simulation if necessary. 

The aluminum antenna is evaporated on top of a sapphire substrate, the thickness of which is 
430 /im for samples 1, 2, 4 and 5 and 500 yum for samples 3 and 6. The contraction of aluminum 
with decreasing temperature leading to a shrinkage of the cavity of about 0.5% and the reduction 
of the permittivity of sapphire by less than a per cent are taken into account [30, 31]. 

The imaginary and real parts of the resulting admittance Y are shown in Fig. 4 for Cj = 0.34ff 
over a range of frequencies spanning three modes. The lowest mode with the largest slope is 
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FIG. 3. (Color online) HFSS model of a 3D-transmon. (a) The 3D resonator with input and output ports. 
These are terminated by 50 Ohm ports, (b) Transparent view of the cavity showing the sapphire substrate. 
Because the electric field is concentrated in the dielectric, a finer mesh is used, (c) and (d) Zoom-ins on 
the antenna placed on top of the substrate. The mesh is finest around the antenna. 




FIG. 4. (Color online) Real and imaginary parts of the admittance Y{u]) = jtoCj — -^jj^ + Yc{ijj) as obtained 
from the HFSS simulation. 
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identified with the qubit mode and the remaining ones with cavity modes, although it must be 
kept in mind that the states corresponding to these modes are ah superpositions of the bare modes. 
With this input, the corrections due to the nonhnearity of the junction are computed as explained 
in the main text. The results of the fit in Cj are given in Table I of the main text and plotted in 
Fig. 5. Details on the measurement of the spectrum can be found in Paik et al. [3]. 
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FIG. 5. (Color online) Low-energy spectrum of six 3D-transmons. Tfieory values (open circles) are obtained 
by fitting Cj for each data set (stars). The green symbols correspond to the — >■ 1 qubit transition, the 
red symbols to the lowest cavity resonance and the cyan symbols to the — >■ 2 qubit transition. The left 
inset show the sub % level relative errors between theory and experiment and the right inset shows the 
fitted values of Cj. 



BLACK-BOX QUANTIZATION WITH MULTIPLE JUNCTIONS 

For simplicity we focus on the dissipationless case. We consider a system with N Josephson 

is) . . (s) 

junctions with bare Josephson energies Ej and charging energies , s ~ 1,...,N, in par- 
allel with a common linear dissipationless but otherwise arbitrary electromagnetic resonator as 
depicted in Fig. 6 (a). The unbiased isolated junctions alone are described by the Hamiltonian 
Hj ~ X^fLi (^^E^\ns)^ — Ej^"^ cos{^ps)j , where Us is the Cooper-pair number operator of the 

s-th junction conjugate to the phase degree of freedom ips, i-C [(/?s,ns] = ih.['S2] A correspond- 
ing iV-port linear circuit, shown in Fig. 6 (b), is then defined by associating a port with each 
junction and replacing the latter with a parallel lumped element LC oscillator with inductance 
l'j^ — {(jy^^Y / e[P and capacitance C'f'' = e^/(2£'p'). Here and in the following 0o = h/{2e) is 
the reduced flux quantum. This corresponds to expanding the cosines in Hj to second order in 
ips- We next consider this linearized circuit classically. 

A quantity of central importance in the following is the A'^-port impedance matrix Z with ele- 
ments Zagi (ijj) ~ Vs [uj] /Is' (w) . Let US choose arbitrarily one reference port k among the N 

ports. By virtue of Foster's theorem [20] Zkki^^) is a purely imaginary meromorphic function and 
can be synthesized by the equivalent circuit of parallel LC oscillators in series shown in Fig. 6 (c). 
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FIG. 6. (Color online) (a) Schematics of A'' JJs (gray (red) boxed crosses) coupled to an arbitrary linear 
circuit (striped disk), (b) Corresponding linearized A'^-port circuit with JJs replaced by parallel LC oscil- 
lators, (c) Foster-equivalent circuit of the impedance Zkk{ijj) of the linearized circuit shown in (b). The 
reference port fc £ {1, . . . , A'^} may be chosen arbitrarily. 



Explicitly 

z,,(.) = f;Lcr + -^) , (12) 

p=i V J^L'p J 

where M is the number of modes and we have adopted the electrical engineering convention of 
writing the imaginary unit as j = —i. This equivalent circuit mapping corresponds, in electrical 
engineering language, to diagonalizing the linearized system of coupled harmonic oscillators. Ac- 
cordingly, the eigen- frequencies cUp = [L^p'^C^^)"^ are determined by the poles of Z^k or more 
conveniently by the real roots of the admittance defined as Yk = Z~^^ and the effective capacitances 

are determined by the frequency derivative on resonance of the latter as = (l/2)ImY^!(wp). 
Note that [20] ImY^'(w) > 0. The Lagrangian of the system can be written as 

p=i \ ) 

in terms of the normal (flux) coordinates (j)^\t) = f^e^'^^* + (/p associated with each of 
the equivalent LC oscillators. From this, we can immediately write the Hamiltonian function of the 
equivalent circuit as Ho = 2 X]^i(/p )*(^p'^^)~^/p ' where the subscript indicates that we consider 
the linear circuit (Fig. 6 (b)). Note that the eigen-frequencies do not depend on the choice of port, 
while the eigenmodes do. Kirchhoff's voltage law implies that up to an arbitrary constant, ^pk{t) = 
4>Q^ J2pLi 4'p'\'t)i where according to Josephson's second relation, ipk{t) = ^ Vfc(r)dr is the 
phase variable of the fc-th (reference) junction with voltage 14. Importantly this simple relation 
holds only for the junction at the reference port k. In order to find the corresponding expressions 
for the other junctions (s ^ k), we notice that the AC voltage amplitude Vs{uj) = jLO(f>i''\uj) at 
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frequency w generated across port s in response to a current with amplitude Is' (w) applied at port 
s' is given by T4(a;) = Z^s' (w)/^' (w). Hence we have ipi''\uj) = {Zsk{'^) / Zkk{^))^k{<-^)- Combining 
this with the above we find that 

^'^\t) = g (/^^e^--.* + (4^)*e-^--^*) . (13) 

Quantization is achieve in the canonical way [23, 24] by replacing the flux amplitudes of the 
equivalent oscillators by operators as 



with the dimensionless bosonic annihilation (creation) operators ap (ap. Direct substitution yields 

the Hamiltonian Hq = huiiajai of AI uncoupled harmonic oscillators (omitting the zero point 
energies) and the Schrodinger operator of phase of the l-th junction is 



M 



(fe) 



ing that (jTs ' = (j)o(ps- 

the dependence on the reference port. Accordingly the root mean square fluctuation of the 
flux of junction s in the multi-mode Fock state 712, . . . , jt-m) is given by \J {{<j)s^^)'^) = 



This is Eq. (7) of the main text using that c^i'^^ = 4'o(ps- The superscript makes explicit 



The anharmonic terms generated by the non-linearity of the Josephson inductance, necessary to 
build a qubit, are included by expressing the higher order terms in the expansion of the cosine in 
the harmonic basis. Including up to the quartic terms we obtain explicitly after normal ordering 

H = Hq - ^ -fpp' (2alap' + a^aj,, + flpCp'^ (16) 
pp' 

N 

- i^PP'm' (Saji^/aqag' + 4ataj^,a;^ag' + AalaprngGg' + apaprngUg' +44/aJaJ,) +^0{ip'^J , 

pp' qq' s—1 

with coefficients ^pp'gg' = J27=i rjTTTyfspCsp'Csg^sg' and jpp' = G^^^i/^ggpp' where, choosing the 



first port as the reference port, — zl\[uj''] \J ^ip ■ Treating the (p"^ nonlinearity in first order 
perturbation theory, one obtains the expressions for the energy, generalized chi-shift and generalized 
anharmonicity given by Eq. (9)) of the main text. 
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